Question 2

if (!require("pacman")) 
  install.packages("pacman")

pacman::p_load(tidyverse,
               dplyr,
               ggplot2,
               tidyr,
               scales,
               readr,
               ggthemes,
               rnaturalearth,
               sf,
               leaflet
               )

# Loading the main datasets for visualisations
livestock_data <- read.csv("./data/Production_Crops_Livestock_E_All_Data/Production_Crops_Livestock_E_All_Data.csv")
area_codes <- read.csv("./data/Production_Crops_Livestock_E_All_Data/Production_Crops_Livestock_E_AreaCodes.csv")
elements <- read.csv("./data/Production_Crops_Livestock_E_All_Data/Production_Crops_Livestock_E_Elements.csv")
item_codes <- read.csv("./data/Production_Crops_Livestock_E_All_Data/Production_Crops_Livestock_E_ItemCodes.csv")

# Identifying livestock commodities that we want
livestock_items <- c("Beef and Buffalo Meat, primary", "Fat of pigs", "Meat, Poultry", "Sheep and Goat Meat", "Milk, Total")

# Filtering for livestock production and trade data
livestock_data_filtered <- livestock_data %>%
  filter(Item %in% livestock_items,
         Element %in% c("Production", "Area harvested", "Yield")) %>%
  select(Area, Item, Element, starts_with("Y"))

columns_to_pivot <- grep("^Y[0-9]{4}$", names(livestock_data_filtered), value = TRUE)

# Reshaping the data to long format
livestock_data_long <- livestock_data_filtered %>%
  pivot_longer(
    cols = all_of(columns_to_pivot),
    names_to = "Year",
    values_to = "Value"
  ) %>%
  mutate(Year = as.numeric(gsub("Y", "", Year))) # To convert Year to numeric

# Remove rows with missing Year values if there are any
livestock_data_long <- livestock_data_long %>%
  filter(!is.na(Year))

# Aggregate data by year and element
aggregated_data <- livestock_data_long %>%
  group_by(Year, Element) %>%
  summarise(Total_Value = sum(Value, na.rm = TRUE), .groups = "drop")

# Step 4: Decode area names and item names
area_decoded <- area_codes %>%
  rename(Area_Code = `Area.Code`, Area_Name = Area)

item_decoded <- item_codes %>%
  rename(Item_Code = `Item.Code`, Item_Name = `Item`)

# For plot 2
commodity_trends <- livestock_data_long %>%
  filter(Element == "Production") %>%
  group_by(Year, Item) %>%
  summarise(Total_Production = sum(Value, na.rm = TRUE), .groups = "drop")

# For plot 3
# Load Natural Earth data for country boundaries
world <- ne_countries(scale = "medium", returnclass = "sf")

# Prepare data for map-based visualizations
# Summarize livestock production by country
map_data <- livestock_data_long %>%
  filter(Element %in% c("Production", "Export Quantity", "Import Quantity")) %>%
  group_by(Area, Element) %>%
  summarize(Total_Value = sum(Value, na.rm = TRUE)) %>%
  ungroup()

# Merge map data with country boundaries
map_data <- map_data %>%
  left_join(world, by = c("Area" = "admin"))

# Create separate datasets for production, export, and import
map_production <- map_data %>%
  filter(Element == "Production")

map_trade <- map_data %>%
  filter(Element %in% c("Export Quantity", "Import Quantity")) %>%
  mutate(Trade_Type = ifelse(Element == "Export Quantity", "Export", "Import"))

#For Interactive plot 4
# Ensure `geometry` column is in correct SF format
map_production_sf <- map_production %>%
  st_as_sf()

population_data <- data.frame(
  Area = unique(map_production$Area),
  Population = runif(length(unique(map_production$Area)), 1e6, 1e8) # Simulated population
)

# Merge population data to calculate per capita production
map_production <- map_production %>%
  left_join(population_data, by = "Area") %>%
  mutate(Per_Capita_Production = Total_Value / Population)

# Convert to sf object
map_production_sf <- map_production %>%
  st_as_sf()

# Step 2: Define color palette
production_pal <- colorNumeric("YlGn", domain = map_production_sf$Total_Value)
per_capita_pal <- colorNumeric("Blues", domain = map_production_sf$Per_Capita_Production)

# Load required libraries
library(ggplot2)
library(dplyr)

# Filter the data to only include "Production" for the Element
aggregated_data_production <- aggregated_data %>%
  filter(Element == "Production")

# Create the line plot for production only
ggplot(data = aggregated_data_production, aes(x = Year, y = Total_Value / 1e6, color = Element)) +
  geom_line(linewidth = 1.2) +
  labs(
    title = "Global Trends in Livestock Production",
    x = "Year",
    y = "Volume (Millions of Tons)",
    color = "Element"
  ) +
  theme(
    plot.title = element_text(size = 22, face = "bold", hjust = 0.5),
    axis.title.x = element_text(size = 19),
    axis.title.y = element_text(size = 19),
    axis.text = element_text(size = 12),
    plot.caption = element_text(size = 14, hjust = 1, color = "grey30"),
    legend.position = "right",  # Legend on the right
    legend.direction = "vertical",  # Stack legend vertically
    legend.title = element_text(size = 19),
    legend.text = element_text(size = 18),
    panel.background = element_rect(fill = "white", color = NA),  # Set white background
    plot.background = element_rect(fill = "white", color = NA),  # Ensure white plot background
    panel.grid.major = element_line(color = "grey90"),  # Subtle grid lines
    panel.grid.minor = element_blank()  # Remove minor grid lines
  ) +
  scale_y_continuous(labels = scales::comma)

# Plot 2: Trends in individual livestock commodities (Production Only)
ggplot(data = commodity_trends, aes(x = Year, y = Total_Production / 1e6, fill = Item)) +
  geom_area(position = "stack", alpha = 0.7) +
  labs(
    title = "Trends in Production of Livestock Commodities",
    x = "Year",
    y = "Production Volume ( Millions - Tons )",
    fill = "Commodity"
  ) +
   theme(
    plot.title = element_text(size = 22, face = "bold", hjust = 0.5),
    axis.title.x = element_text(size = 19),
    axis.title.y = element_text(size = 19),
    axis.text = element_text(size = 12),
    plot.caption = element_text(size = 14, hjust = 1, color = "grey30"),
    legend.position = "right",  # Legend on the right
    legend.direction = "vertical",  # Stack legend vertically
    legend.title = element_text(size = 19),
    legend.text = element_text(size = 18),
    panel.background = element_rect(fill = "white", color = NA),  # Set white background
    plot.background = element_rect(fill = "white", color = NA),  # Ensure white plot background
    panel.grid.major = element_line(color = "grey90"),  # Subtle grid lines
    panel.grid.minor = element_blank()  # Remove minor grid lines
  )

# Plot 3: Livestock Production on a World Map
ggplot() +
  geom_sf(data = world) +
  geom_sf(data = map_production, aes(geometry = geometry, fill = Total_Value), color = "black") +
  scale_fill_viridis_c(name = "Production (tons)", option = "C") +
  scale_fill_gradient(
    low = "khaki1", 
    high = "green4"
  )+
  labs(
    title = "Global Livestock Production",
    subtitle = "Production volumes (tons) by country",
    caption = "Source: FAO Livestock Data"
  ) +
  theme_minimal()

# Enhanced Leaflet Map
leaflet() %>%
  # Add a base map
  addProviderTiles(providers$CartoDB.Positron) %>%
  
  # Add Total Production Layer
  addPolygons(
    data = map_production_sf,
    fillColor = ~production_pal(Total_Value),
    color = "black",
    weight = 1,
    fillOpacity = ~pmax(Total_Value / max(map_production_sf$Total_Value, na.rm = TRUE), 0.2),  # Dynamic opacity
    label = ~paste(
      "<strong>Country:</strong> ", Area, "<br>",
      "<strong>Total Production:</strong> ", scales::comma(Total_Value), " tons"
    ),
    popup = ~paste(
      "<strong>Country:</strong> ", Area, "<br>",
      "<strong>Total Production (tons):</strong> ", scales::comma(Total_Value), "<br>",
      "<strong>Population:</strong> ", scales::comma(Population), "<br>",
      "<strong>Per Capita Production (tons):</strong> ", round(Per_Capita_Production, 2)
    ),
    highlightOptions = highlightOptions(
      weight = 3,
      color = "blue",
      bringToFront = TRUE
    ),
    group = "Total Production"
  ) %>%
  
  # Add Per Capita Production Layer
  addPolygons(
    data = map_production_sf,
    fillColor = ~per_capita_pal(Per_Capita_Production),
    color = "black",
    weight = 1,
    fillOpacity = ~pmax(Per_Capita_Production / max(map_production_sf$Per_Capita_Production, na.rm = TRUE), 0.2),  # Dynamic opacity
    label = ~paste(
      "<strong>Country:</strong> ", Area, "<br>",
      "<strong>Per Capita Production:</strong> ", round(Per_Capita_Production, 2), " tons"
    ),
    popup = ~paste(
      "<strong>Country:</strong> ", Area, "<br>",
      "<strong>Per Capita Production (tons):</strong> ", round(Per_Capita_Production, 2), "<br>",
      "<strong>Total Production (tons):</strong> ", scales::comma(Total_Value), "<br>",
      "<strong>Population:</strong> ", scales::comma(Population)
    ),
    highlightOptions = highlightOptions(
      weight = 3,
      color = "green",
      bringToFront = TRUE
    ),
    group = "Per Capita Production"
  ) %>%
  
  # Add Layer Control
  addLayersControl(
    baseGroups = c("Base Map"),
    overlayGroups = c("Total Production", "Per Capita Production"),
    options = layersControlOptions(collapsed = FALSE)
  ) %>%
  
  # Add Legend for Total Production
  addLegend(
    pal = production_pal,
    values = map_production_sf$Total_Value,
    title = "Total Production<br>(in tons)",
    position = "bottomleft",
    group = "Total Production"
  ) %>%
  
  # Add Legend for Per Capita Production
  addLegend(
    pal = per_capita_pal,
    values = map_production_sf$Per_Capita_Production,
    title = "Per Capita<br>Production (tons)",
    position = "bottomright",
    group = "Per Capita Production"
  )